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Abstract. In this paper we present an optimal control approach modeling fast exit scenar- 
ios in pedestrian crowds. In particular we consider the case of a large human crowd trying to 
C^ exit a room as fast as possible. The motion of every pedestrian is determined by minimizing 

' ' a cost functional, which depends on his/her position, velocity, exit time and the overall den- 

—T" sity of people. This microscopic setup leads in the mean-field limit to a parabolic optimal 

control problem. We discuss the modeling of the macroscopic optimal control approach and 
t^ show how the optimal conditions relate to Hughes model for pedestrian flow. Furthermore 

^^ we provide results on the existence and uniqueness of minimizers and illustrate the behavior 

'^N of the model with various numerical results. 
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1. Introduction 



Mathematical modeling of human crowd motion such as streetway pedestrian flows or the 
evacuation of large buildings is a topic of high practical relevance. The complex behavior 
,^ of human crowds poses a significant challenge in the modeling and a majority of questions 

5^ still remains open. Consequently it receives increasing attention, and in particular in the 

last years strong development in the mathematical literature becomes visible. A variety of 
different approaches has been proposed and partly analyzed, which can roughly be grouped 
as follows: 

^ (1) Microscopic force-based models such as the social force model (cf. |26t [25| I14j ) and 

^^ their continuum limits. 

^ (2) Cellular automata approaches (discrete random walks with exclusion principles p^[6]) 

l/^ and their continuum limits (cf. |10])- 

■^ (3) Fluid- dynamical and related macroscopic models (cf. |15 | l32 [[roi [5]). among which in 

O particular the Hughes model and its variants received strong attention (cf. [^ [T71 [U 

cn [52]). 

. . (4) Microscopic optimal control and gavae approaches (cf. [571 EB]) and mean field game 

^ models (cf. |301ll8j). which can be considered as the appropriate continuum limit (cf. 

X lamiii). 

For evacuation scenarios an additional complication in the basic model arises, namely how the 
boundaries (i.e. doors, impermeable walls) and the goal of the pedestrians to leave the room 
as quickly as possible are modeled. A canonical approach is to use a potential force (or drift 
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in random walk models) based on the distance function to the doors. This means that the 
potential (j) solves some kind of Eikonal equation in the domain Q with zero Dirichlet boundary 
conditions on the part Te C 9il, which corresponds to the door. The appropriate statement 
of the Eikonal equation and its boundary conditions in a crowded situation is a delicate issue, 
which is mainly done in an ad-hoc fashion (even for detailed microscopic models) and has led 
to different directions. Some models simply use the distance function itself, i.e. the viscosity 
solution of 

|V(^| = 1 inQ., (1.1) 

which corresponds to an overall optimization of the evacuation path independent of other 
pedestrians. This might be questionable in crowds due to limited visibility or global changes 
of the evacuation path to avoid jamming regions. Other models use variants of the Eikonal 
equation that incorporate the density p of pedestrians. The celebrated Hughes model, cf. 
p8] , uses an equation of the form 

/(p)|V0| = l inf], (1.2) 

where / is a function introducing saturation effects such as f(p) = pmax — P, where pmax is a 
maximal density. 

In this paper we introduce a rather unifying mean-field game approach to the modeling of 
evacuation scenarios, which links several of the above mentioned models and provides further 
understanding of several issues. This mean field framework provides interpretations of differ- 
ent forms of the Eikonal equation and covers the modeling of 'realistic' boundary conditions. 
We provide a detailed mathematical analysis and numerical simulations of our approach and 
discuss interesting special cases and limits, in particular its connection to Hughes' model. 
This paper is organized as follows: In Section [2] we present the modeling on the microscopic 
level and its mean-field limit. Furthermore we discuss its relation to Hughes' model and 
the appropriate choice of boundary conditions. Section [3] is devoted to the analysis of the 
macroscopic optimal control approach. Finally we present a steepest descent approach to 
solve the parabolic optimal control problem in Section |4] and illustrate the behavior with 
several numerical examples. 

2. Mathematical Modeling 

In the following we discuss a model paradigm based on the idea of fast exit. At the single 
particle level, this yields a classical or stochastic optimal control problem, which we refor- 
mulate as PDE-constrained optimization model for the particle density. From the optimality 
conditions we then obtain a version of the Eikonal equation as the adjoint problem. 

2.1. Fast Exit of Particles. Let us start with a deterministic particle (of unit mass) trying 
to leave the domain il. as fast as possible. Let X = X{t) denote the particle trajectory and 

Te^itiX) = sup{t > I X{t) E n}. 

Then it makes sense to look for a weighted minimization of the exit time T^xit and the kinetic 
energy, i.e. 

/•Te: 

subject to X(t) = V{t), X{0) = Xq, where a > encodes the weighting of the fast exit. 



^ ^ ''" |y(t)|2 dt + ^TeMX) ^ mm, (2.1) 
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Introducing the Dirac measure ii = ^xit)^ ^^^ ^ final time T sufficiently large, we can rewrite 
the functional (2.1) as 

rT 



/t(^,w) 



Jn 



\v{x,t)\ dfi dt + 



a 



Jn 



dfj, dt, 



subject to 

dtH + V- (/iv) = 

and initial value iJ,\t=o = Sxo- Note that for T > T^xu 

TexityX) 



(2.2) 
(2.3) 



Jo. 



d6x{t) dt 



and 



Jn 



\v{x,t)\'^dii dt 



Jn 



\v{x,t)\^ dSx(t) dt 



|u(X(t),t)|2 dt, 



which yields the equivalence between the particle formulation (2.1 ) and the continuum version 
(2.2) by the standard Lagrange-Euler transform of the velocity V{t) = v{X{t),t), cf. [271 [36]. 
Next we consider a stochastic particle, which moves according to the Langevin equation. 
Then we obtain 

dX{t) = V{t) dt + a dW{t), (2.4a) 

where W a Wiener process and a the diffusivity. It is natural to consider the stochastic 
optimal control problem 

'1 



E 



Xo 
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. , \V{t)\'dt^-T,xiAX) 
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(2.4b) 



with the random variable X determined by ( |2.4a ) with initial value X^. Reformulating ( |2.4[ ) 
for the distribution already reveals the structure of a mean-field game for the particle density. 
That is, writing dfi = p dx we obtain the minimization functional 



subject to 



It{p,v) = 7: / p{x,t) \v{x,t)\ dx dt + - / p{x,t) dx dt, (2.5a) 

^ Jo Jn 2 Jq J^ 



dtp + V- {pv) = ^Ap, p{x, 0) = pq{x). (2.5b) 



2.2. Optimality, Eikonal Equations, and Mean-Field Games. In the following we pro- 
vide an explanatory derivation of the optimality conditions for the optimization problem 



It{p,v) — )■ min . 

{p,v) satisfying l|2.5b^ 



(2.6) 



We shall derive the optimality condition at a formal level, without providing rigorous details 
on the functional setting and on the boundary conditions. We shall be more precise later on 
a more general model, see subsection |2.6[ 
We define the Lagrangian with dual variable (p = (l){x,t) as 

(2.7) 



Lt(p, V, 4>) = It{p, v)+ I / (dtp + V • {pv) - —Ap)(p dx dt. 

Jo Jn ^ 

For the optimal solution we have, in addition to (2.5b), the following equations, i.e. 

= 9„Lt(p, V, (p) = pv - pV(p, 



(2.^ 
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and 

= dpLrip, V, <P) = ^|z;p + ^ - dt(^ - v ■V(t> - y A</>, (2.9) 

with the additional terminal condition (j){x,T) = 0. Inserting v = \/(j) we obtain the following 
system, which has the structure of a mean field game, i.e. 

2 

^ip + V • (pV^) - — Ap = (2.10a) 

dt<t>+\\Vct^\^+'^/^<P = ^. (2.10b) 

One observes the automatic emergence of a transient viscous Eikonal equation, which is to 
be solved backward in time with terminal condition. The connection to models using the 
distance function as a potential is quite apparent for a time interval [0, S] with S « T and 
(7 = 0. Noticing that the Hamilton-Jacobi equation is solved backwards in time and that the 
backward time is large for t < S, we see that the solution is mainly determined by the 
large-time asymptotics (cf. [29j) solving 

\V4>\^ = c, (2.11) 

for some constant c. Hence the potential becomes approximately a multiple of the distance 
function. 

2.3. Mean Field Games and Crowding. We now turn our attention to a more general 
mean field game, respectively its optimal control formulation. We generalize the terms de- 
pending on the density for reasons to be explained in detail below to obtain 

It{p,v) = 1 I F{p)\v{x,t)\^dxdt + \ [ E{p) dx dt, (2.12a) 

'^ Jo Jn 2 Jo Jq 



subject to 



a^ 



dtp + V-{G{p)v) = -Ap, (2.12b) 

and a given initial value p{x,0) = po{x). The motivation for those terms is as follows: 

• The function G = G{p) corresponds to nonlinear mobilities, which are frequently used 
in crowding models, e.g. G{p) = p(pmax — p) derived from microscopic lattice exclusion 
processes (cf. [1Q\ for a pedestrian case). Nonlinear mobilities have been derived in 
different settings and used also in different applications of crowded motion, e.g. ion 
channels [9l [11] or cell biology |35l [33l [19]. A typical form of G is close to linear 
and increasing for small densities, while saturating and possibly decreasing to zero for 
larger densities. 

• The function F = F{p) corresponds to transport costs created by large densities. 
In particular one may think of F tending to infinity as p tends towards a maximal 
density. As we shall see below however, changes in F can also be related equivalently 
to changes in the mobility. 

• A nonlinear function E = E{p) may model active avoidance of jams in the exit 
strategy, in particular by penalizing large density regions. 



The optimality conditions of (2.12 ) can be (again formaly) derived via the Lagrange functional 

Lt{p, V, 0) = hip, v)+ [ [ (dtp + V- {G{p)v) - ^ Ap)0 dx dt. (2.13) 

Jo Jn ^ 
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For the optimal solution we have 

= 9„Lt(p, V, <P) = F{p)v - G{p)Vcl), (2.14) 

and 

= dpLrip, V, </.) = ^F'{p)\v\^ + '^E'ip) - dt^ - G'{p)v ■ V0 - ^A4>, (2.15) 

with terminal condition (p{x, T) = 0. Inserting v = j^Vcj) we obtain the optimality system 

dtp + V- (^V<A) - ^Ap = (2.16a) 

F{p) 2 

5t0 + 2 (2^ - -^)\^<^? + y Ac/> = -^E\p), (2.16b) 

again a mean field game. 

2.4. Momentum Formulation. In order to reduce ambiguities in the definition of a velocity 
and also to simplify the constraint PDE, we use a momentum or flux-based formulation in 
the following, cf. [2], [5]. The flux density (or momentum) is given by j = G{p)v, hence the 
functional do be minimized in (2.12) becomes 

MpJ) = 1 [ [ TTTT^^ ^* + ^ / / ^(/') ^^ ^*' (2.17a) 

2 Jo Jn H{p) 2 Jq J^ 

with H := ^^, subject to 

2 

dtp + V-j = ^Ap. (2.17b) 



G- 



We see the redundancy of F and G in the above formulation, effectively only H .— p 
determines different cases. We mention that in order to obtain a rigorous formulation, we 



:|2 



replace j^J-^ in (|2.17a|) by 



j! 



K{3, p):=\ if j = 0, H{p) = (2-18) 

+00 if j / 0, H{p) = 0. 

The Lagrange functional for this problem is given by 

Lt{p,v,4>)=It{p,J)+ I [ idtp + V-j-^Ap)cl)dxdt, (2.19) 

Jo Jn ^ 

and the optimality conditions are 



2 

9tp + V • {H{p)V(l)) - y Ap = (2.20a) 

dt<P + ^|V<A|2 + ^A<p = ^E'ip). (2.20b) 

We recall that the rigorous formulation of the problem in terms of the boundary conditions 



will be performed in subsection 2.6 
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2.5. Relation to the Hughes Model. For vanishing viscosity a = the optimahty system 



(2.16) respectively (2.20) has a similar structure as Hughes model for pedestrian flow, see 



[28], which reads as 

dtp-V-{pf{p)V(l)) = (2.21a) 

|V<^| = -^. (2.21b) 

Hughes proposed that pedestrians seek the fastest way to the exit, but at the same time try to 
avoid congested areas. Let />max denote the maximum density, then the function f{p) models 
how pedestrians change their direction and velocity due to the overall density. A common 
choice for example is f{p) = pmax — P- 
To obtain the connection with the Hughes model we choose H(p) = pf(p)'^, E(p) = ap, then 



(2.20) reads as 

dtp + V- {p f{pfV4>) = (2.22a) 

9*0 + ^(/(p) + 2p/'(p))|V<Ap = |. (2.22b) 

For large T we expect equilibration of (j) backward in time. If we further consider cr = 0, then 
for time t of order one the limiting model (very formal) becomes 

dtp + V-{pf{pfV(t>) = (2.23a) 

(/(;.) + 2p/'(p))|V0p = -^, (2.23b) 



which is almost equivalent to Hughes model (2.21[), if we set c = 1. Note that the difference 



is in the term 2f'{p), which is however severe. For f{p) = Pmax — p vie have 

f{p) + 2p/'(p) = Pmax - 3p, 

thus for small densities the behavior is similar, but the singular point is p = ^^y^. A (so far 
partial) attempt to a rigorous mathematical theory for the original Hughes model (|2.21) has 



been provided in [T]. Previous results |17| I15| considered smoothed version of that model. 
El-Khatib et al. first considered a variant with a cost in the eikonal equation not related to 



/(/)), cf. |21| . and finally Gomes and Saude provided a priori estimates for (2.22) in 



2.6. Boundary Conditions. We finally turn to the modeling of the boundary conditions. 



which we have neglected in the computations above. To do so we consider (2.17) on a bounded 
domain Q C M*^, d = 1,2, corresponding to a room with one or several exits. We assume 
that the boundary dQ is split into a Neumann part Ftv C Ofi, modeling walls, and the exits 
r_E C dCl, do, = Tjy U r^; and F^r n F^; = 0. We now formulate natural boundary conditions 
on the density, which will lead to adjoint boundary conditions for (p. 
On the Neumann boundary F^r we clearly have no outflux, hence naturally 

2 

(-YVp + j)-n = 0, 

where n denotes the unit outer normal vector. The exit part is more difficult. Here the outfiux 
depends on how fast people can leave the room. If we denote the rate of passing the exit by 
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j3, then we have the outflow proportional to j3p. Hence, we arrive at the Robin boundary 
condition 

,2 



a 



{-—Vp + j)-n = ^p. 



(2.24) 



Then boundary conditions for the adjoint variable (p can be calculated again from the opti- 
mality conditions. If we define the Lagrangian in this case we obtain: 

LT{p,j,(t)) = It{p,J) + {dtP + V-j- —Ap)^ dx dt 

Jo Jn ^ 

= hip, j)+ f f p{-dt(t) - ^ A0 - j ■ V(/>) dx dt 
Jo Jn ^ 

rT [■ 2 2 

+ / ( Vp • ncj) H pV4> ■ n + j ■ ncj)) ds dt 

Jo JdVE 2 2 

-pSJ(j) ■ n ds dt 



+ 



T . ^2 



JdVM 2 



It{p,J) + 

Jo 



p{-dt(f) 



n 



T 



a 



A0 - j • V(/)) dx dt 

T 



+ / / (/3/90 H pV<j) ■ n)ds dt + 

'o JdVE 2 Jo JdTr 



a 



p^cf) ■ n ds dt. 



(2.25) 



In this form we see that the optimality condition with respect to p, i.e. dpLx = 0, yields the 
following boundary conditions for the adjoint variable <j) in addition to the adjoint PDE: 



2 2 

— V(/> • n + /?</> = on r^; and — V0 • n = on Tjy. 



(2.26) 



One observes that our model provides meaningful and easily interpretative boundary condi- 
tions for the density p as well as for the adjoint variable (j)- A homogeneous Dirichlet boundary 
condition for (p on doors arises only in the asymptotic limit /3 — )• c». In this case one also 
needs to specify a homogeneous Dirichlet boundary condition for p for consistence. This is in 
contrast to previous models, which frequently used Dirichlet boundary conditions for 0, but 
arbitrary Dirichlet values for p or an outflux boundary condition as above with finite /3. 



In the case of Hughes model (2.21 ) commonly used boundary conditions are chosen as follows. 
The potential (j) is assumed to satisfy homogeneous Dirichlet boundary conditions, i.e. <J3 = 0, 
while the density p either satisfies a Dirichlet boundary conditions of type p = po oi an 



of (2.21) and (2.22): 



outflux condition as in (2.24). There are two possibilities to compare the behavior of solutions 



= pd and homogeneous Dirichlet 
This, however, is not realistic in 



Either by prescribing a density at the exit, i.e. p 

boundary conditions for the adjoint, i.e. = 0. 

practical applications. 

Or by choosing outflux boundary conditions for the density p, as in (2.24) with a 



large proportionality constant /3, and the corresponding no-flux boundary conditions 
(2.26) in the for (p. The latter choice for the potential (p is not compatible with (2.21), 



but approximates the homogeneous Dirichlet boundary condition (typically chosen in 
Hughes model) for /3 large. 
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Note that no flux boundary conditions for tlie density p, i.e. /3 = in (2.24) would result in 



a homogeneous Neumann boundary conditions (2.26) for (p. Such boundary conditions are 



not compatible with the Eikonal equation (2.21b) in Hughes model and contradict the initial 
modeling assumptions. 

3. Analysis of the Optimal Control Model 
3.1. Existence of Minimizers. In this section we discuss the existence and uniqueness of 



minimizers for the optimization problem (2.12), respectively the reformulation (2.17). In 



order to avoid any ambiguity with the formulations we can simply set 

F = G = H, 

which will remain as a standing assumption throughout the rest of this section. Let i7 be a 
bounded domain in M , d = 1, 2,. We make the following basic assumption on F and E: 

(Al) F = F{p) G C^{^), F bounded, E = E{p) G C^{W) and F{p) > 0, E{p) > for 

per. 

Existence of minimizers is guaranteed if 
(A2) E = E{p) is convex. 



As recalled in subsection 2.5, models with maximal density are of great importance in the 
applied setting, therefore we shall consider the special class of mobilities F which yields a 
density uniformly bounded by a maximal density. Let us fix the maximal density to be 
Pmax > 0. Let us denote by T = [0,pniax]- Typically pmax = 1- In order to force the 
minimizers of the problem to satisfy p G T, we extend F to F{p) = on p G M \ T. This 
defines the next assumption on F, i.e. 

(A3) F(0) > if p G T and F = otherwise. 

Note that assumption |( A3)"] is satisfied in particular by F{p) = pf{p) in the Hughes model in 
subsection |2.5[ To show uniqueness we need another additional assumption for F, namely 

(A4) F = F{p) is concave. 

We consider the optimization problem on the set V x Q, i.e. It{p-, t') : ^ x Q — t- M, where V 
and Q are defined as follows 

V = L^{id,T;H^{n))r\H^{Q,T;H-^{n)) andQ = L'^{Qx (0,r)). (3.1) 

Hence the optimization problem (|2.12[) reads as 



a^ 



min It{p,v) such that dtp = — Ap — V • {F{p)v), (3.2) 

(p,v)£VxQ 2 

respectively in momentum formulation 

2 

min It{p,j) such that dtp = — Ap — V • (j). (3-3) 

ip,j)eVxQ 2 

In order to have the differential constraint well defined in the prescribed functional setting, 
and to incorporate the Robin boundary conditions in the statement of the problem, we shall 
provide a more rigorous definition of the minimization problem below. 

For the existence proof in the case of non-concave F, we introduce another formulation based 
on the rather nonphysical variable 

w = y^F{p)v, 
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which has been used aheady in the case of hnear F in [8j. Then the functional It can be 
re-written as 

1 



J{p,w) 



2 Jo Jn 
and the optimization problem formally becomes 



{\w\^ + E{p)) dx dt, 



a 



min J(/9, w) such that dtp = — Ap — V • {y F{p)w). 
{p,j)eVxQ 2 



(3.4) 



(3.5) 



In order to make the relation among the problems (3.2), (3.3), and ( |3.5[ ) rigorous, we need to 
extend the domain of the velocity v to 



Qp := {v measurable | \/F{p)v G Q}. 
Moreover, for given p we define an extension mapping w £ Q to v £ Qp via 



w{x) 



Rp{w){x) :-- 



else. 



iiFipix))^0 



(3.6) 



(3.7) 



With this notation, we easily get (we omit the details) 

Lemma 3.1. Let | (A 1)\ be satisfied and let a > 0. Then the following relations hold: 

• Ifip,v) G VxQp solves ( fO] ), then {p, ^jF{p)v) G FxQ solves ( [3^ . If{p,w) G V y.Q 
solves (3.5), then {p, Rp{w)) £VxQp solves (3.2). 

• // {p,w) e V X Qp solves ( [3l5| ), then {p, ^jF{p)w) ^VxQ solves ( [3^ . // {p,j) G 
V xQ solves (3.3), then {p,Rp{j)) €V x Qp solves (3.5). 

As a consequence of the above Lemma, we see that if one of the models has a minimizer with 
F(p) uniformly bounded away from zero, it is in 1/ x Q and minimizes all three models. We 
are now ready to state the minimization problem (3.2) and its reformulation (3.5) rigorously, 
in a way to incorporate the Robin boundary conditions, and to take into account the proper 
functional setting. 

Definition 3.2 (Minimisation problem, 1st formulation). Let po £ L'^{^). A pair {p,v) G 
V X Qp is a weak solution to the minimization problem (3.2) with initial condition po, if 
p(0) = Po and 

{dtp, ^P)H-Km + / {^^P - Hp)v) ■Vil^dx = - I pp^ ds, (3.8) 

Jn ^ Jte 

for all ij G H^n), and if 

It{p,v) = min {lTip,v), : {p,v) £V x Q, (p,-y) satisfy ([3^} . 



A similar re- formulation of (3.3) in terms of the momentum variable j replacing v can be also 
deduced. 

Definition 3.3 (Minimisation problem, 2nd formulation). Let po £ L'^{0,). A pair (p, j) G 
V X Qp is a weak solution to the minimization problem (3.3) with initial condition po, if 
p(0) = po and 

{dtp, ^)H~\m + I (^Vp -j)-V^dx = - f Ppi; ds, (3.9) 

Jn ^ Jfe 



MARTIN BURGERi, MARCO DI FRANCESCO^, PETER A. MARKOWICH^, AND MARIE-THERESE WOLFRAM* 

for all V G H\n), and if 

It{pJ) = min|7T(p,j), : (pJ) &V xQ, {p,j) satisfy ([3^| , 
where 

It{p,J) 

and K is defined as in ( |2.18 ). 
Similarly, on the variables setting {p, w) we get: 

Definition 3.4 (Minimisation problem, 3rd formulation). Let po S L'^{Q). A pair (p^w) £ 
y X Q is a weak solution to the minimization problem (3.5) with initial condition pQ, if 
p{0) = po and 

^2 



\ I I K{j, p)dx dt + \ I I E{p) dx dt, 
■^ Jo Jn ^ Jo Jn 



,cr 



{dtp, i')H-\H^ + / (y Vp - ^/F{p)w) • VV^ dx 



Ppip ds, 



(3.10) 



for all V G H^i^), and if 

Jxipjw) = mini Jt{p,w), : {p,w)£VxQ, {p,w) satisfy (3.10)}. 

We now prove our main existence result in the formulation using w, which allows also for 
non-concave F, and relies on a positive viscosity a. First, we provide an a-priori estimate: 



Lemma 3.5. Let pQ e L^(i7). Let (Al) and (A2) be satisfied and let a > 0, (5 > 0. Let 
w ^ Q and let p ^V he a weak solution of 



,(y 



{dtp, iP)H-i,m + J ( yVp - y^F{p)w) • V^ dx 



(3 pip ds, 



(3.11) 



for all ip G H^{Q,). Then there exist constants Ci,C2 > depending on F, a, Q and T only, 
such that 

WpWv < Ci\\w\\q + C2. (3.12) 

Proof We use the test function tp = p{t) and obtain 



Id^ 
2Jtjn 



p dx 



\Vp\'^dx + / ^jF{p)wVp dx- j I3p^ ds. 
2 ./ ./o Jyt. 



Since F = F{p) is bounded, the second to last term is non positive, and w{t) £ L^(i7). This 
yields after integration in time 



2 , <^ 
L2 + 



1 



2 IIpWIIl2 + y /„ WPllmdT < ci j^ \\w\\l2 dr + -WpoWh^n)^ 
with a constant ci € M^ depending only on a and F. For the time derivative we deduce that 

\\dtp\\H-Hn)= sup \{pt,ip)\< sup [c2||p||hi||V'IIhi +c3||w||l2||vV'IIl2]> 

where the first term on the right-hand side includes the gradient terms in the domain as well 
as the boundary term via a trace theorem. Squaring and integrating in time finally yields the 
desired estimate for ll/oHy. D 

The following technical Lemma is needed to ensure that assumption |(A3)] gives p € T. We 
provide the details for the sake of completeness, although the result is quite standard. 
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Lemma 3.6. Assume p and w are as in the assumption of Lem^ma 3.5. and assume further 
that F satisfies l(A3)\ Then, p{-,t) G T for all t £ {0,T] if pq{x) G T. 

Proof. Assume first w G C^, and let Tjsip) be a C^'^ regularization of the positive part 
(p)+ = max{/9, 0} as (5 \ 0, such that rj'Ap), v'sip) ^ Oi ^^'^ v'sip) = jXo<p<<5- We can use the 
test function r]'{p — /Omax) in Definition |3.4| to get 



^ / VsiP - /Jmax) dx = -y / Tj'g {p - 


-Pmax) Vpp dx 


+ / 'JP{p)wr}'^{P - Pmax) -Vp dx- 

Jn 


-/3 / prj'sip-p 


< / y/F{p)wr]s{p - pmax) • V/9 dx 

Jn 




= - Cs{p)diyw dx, 





(3.13) 



'n 



where 



^s{p) ■■= ^/F{p)r]'l{z- p^^y)dz. 

Now, it is clear that ^^(/o) = if p < Pmax- Moreover, for Pmax < P < Pmax + 5 we have 

Up) = \f VW)dz < c, 

since F is uniformly bounded. Finally, in the case p > pmax + S, we have 

]^ fPmax + S 

Csip) = ^ / ^/F{^dz ^ V^(Pmax) = 0, 

as (5 \j 0. Combining all these assumptions, we get ^sip) uniformly bounded and such that 
(,sip) ~^ almost everywhere on fi x [0,r]. Since w G C^, we can integrate in time in (3.13) 
and send 5 \ to get 

{p{x,t) - p^ax) + dx < / {po{x) - pmax) + dx = 0, 

n Jo. 

which gives the assertion. The general case w £ Q can be recovered by a standard approxi- 
mation argument, we omit the details. D 



Theorem 3.7 (Existence in the general case). Let po G L (0). Let (Al) and (A2) be satisfied 
and let o" > 0. Then the variational problem (3.5) has at least a weak solution {p,w) ^VxQ 
with initial condition po ^^ the sense of Definition 3.3. If in addition (A3) is satisfied, then 
PGT. 

Proof. We show the existence of a minimizer by the direct method. Let {pk^Wk) G V x Q 
be a minimizing sequence. Then we can assume without loss of generality that there exists a 
constant c, such that 

/ {{wkl"^ + E{pk)) dx dt < c. 
Jo Jvt 

Clearly, Wk is bounded in L'^{Q. x (0, T)) and hence has a weakly convergent sub-sequence with 

limit w, again denoted by Wk- Now we use Lemma 3^ to obtain boundedness of pk G V , from 

which we can extract another weakly convergent sub-sequence. From the Lemma of Aubin 
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and Lions, cf. |33j, we have a compact embedding of pk € L^(il x (0,T)). Then the continuity 
and boundedness of F impUes that yj F{pk) converges strongly in L?{Q. x (0, T)) and the hmit 
yjF{p) is bounded. Hence, the product \/F{pk)wk converges weakly in L^(Ox (0, T)) and the 
limit yjF{p)w is in L^(0 x (0,T)), which allows to pass to the limit in the weak formulation 
of the constraint equation. Thus, (p, w) is admissible, and the weak lower semicontinuity of 
the objective functional implies that it is indeed a minimizer. Note that the convexity of E 
is crucial in this case, as it provides lower semi continuity of j^ E{pjS)dx. Assume now that 
|(A3)| is satisfied. Then, due to Lemma 3.6 for all k we have pk S T, which can be passed 



to the limit (possibly by extracting a further subsequence converging almost everywhere) to 
obtain the desired result. D 

We now provide an alternative argument, which can be extended to the non- viscous case 
for zero-flux flux boundary conditions, i. e. for /3 = 0. The argument below provides less 
regularity for the minimizer compared to the previous result, but the proof is much shorter. 



Theorem 3.8 (Existence for concave mobility). Let \(A 1)\ \ (A 2)\ \ (A 3)\ and \(A4)\ be satisfied 



and let a > 0. Then the variational problem (3.3) has at least one minimizer {p,j) G L°°(0 x 



(0,r)) X Q such that p{x) G T for almost every x G il. If E is strictly convex the minimizer 
is unique. 

Proof. Let {pk,jk) be a minimizing sequence. Since F{pk) is bounded, jk is uniformly bounded 
in Lp'iQ X (0,T)) and we can extract a weakly convergent sub-sequence with limit j. By the 
boundedness of the objective functional, we have F{pk) > almost everywhere, which gives 
Pk{x,t) G T for all k and for all (x,t). From the latter we can extract a sub sequence, again 
denoted by pk such that pk -^* p in L°°(il x (0, T)), respectively weakly in L^ for p G (1, oo). 
Since the constraint PDE is linear we can easily pass to the limit in a weak formulation and 
see that 

2 
5tP + V • (i) = y Ap, 

holds in a weak sense. Moreover, the convexity of the objective functional guaranteed by 
|(A2)[ |(A3)[ and |(A4)| (recall that l/F{p) is extended to -Foo when F{p) = 0) implies weak 



semi-continuity an hence {f),j) is a minimizer. The uniqueness follows from a standard strict 
convexity argument. D 



Remark 3.9. The result in Theorem \3.^ can he extended to the case a = in case the 
boundary conditions are posed with (3 = 0, i. e. no exits at the boundary. In this case, the 
boundary condition can he trivially passed to the limit k — )• +oo above, whereas we are not 
able to close such argument in case /S > 0. This case requires a more refined analysis based 
on the theory of nonlinear conservation laws with boundary conditions, cf. [1] . 

Remark 3.10. We mention that under the above assumptions one might expect Gamma- 
convergence of the optimal control problems as o" — )• 0. The lower semicontinuity arguments 
can be used to verify the lower bound inequality, however it is so far unclear how the upper 
bound inequality can be verified. 

Note that the convexity constraint | ( A2 ) | can be weakened, i.e. the function E has to be weakly 
lower semicontinuous in suitable function spaces to guarantee existence of minimizers. Since 
we only consider convex functions E, we did not state this more general result. 



MEAN FIELD GAMES WITH NONLINEAR MOBILITIES IN PEDESTRIAN DYNAMICS 



13 



3.2. Adjoint Equations and Optimality System. In the following we investigate the 
existence of adjoints, i.e. the optimality condition with respect to the state variable p: 



Proposition 3.11. Let assumption \ (A 1)\ and \(A2)\ be satisfied and let p be such that H{p) > 7 
for some 7 > 0. Then the adjoint equation of problem (3.3), i.e. 



dt 



(/)(x,r) 



E'{p) 



2IJI p2 



(3.14a) 
(3.14b) 



with boundary conditions (2.26) and terminal data 



e L9(0, T; W^''i{n)) with q < ^^^^ 



at t = T has a unique solution 



N+l- 



Note that the right hand side of (3.14) is only in L^{U x (0, T)). Existence and uniqueness 
of solutions (p follows from the results of Boccardo and Gallouet, cf. [7j. However this result 
does not provide the necessary regularity for (j) to define the Lagrange functional Lx properly. 
Finally we investigate the full optimality system, where we eliminate j = F{p)V(j) and obtain 



dtp 

dt(t> 



a 



^ Ap - V • (F(p)V0) 
-^A<p+^-E'{p)-^-F'{p)\VcP\\ 



with boundary conditions 

,2 



a 



Vp + F{p)V(t)) -n 



on Fat 
Pp on F^ 



a 



-'V(p ■ n 




(3.15a) 
(3.15b) 

(3.15c) 



To show uniqueness of the forward backward system (3.15) we follow an approach proposed 
by Lasry and Lions for mean field games, cf. |31j . 



Theorem 3.12 (Uniqueness for the optimality system). For a fixed initial condition po £ 
L^iyi), there exists a unique weak solution 

{p,(t)) G L\0,T;H\n)) x L\0,T;H\n)) 



to the optimality system (3.15) with boundary conditions (3.15c). 



Proof Assume there exists two classical solutions {(j)i,pi) and {(j)2,P2) to the optimality sys- 
tem (2.16). We denote the difference between the solutions hy p = pi — p2 and (p = cpi — (p2, 



which satisfy 
dtp 
dt^ 



a 



Ap-V- (F(pi)V0i - F{p2)V(j)2) 



.^A^+^-E'{p^)-^-E\p2) 



V(/,i)|V0i|2 + 1f'(p2)|V<A2| 



Note that the differences p and (j) satisfy the following boundary conditions on Fg: 

0. 



2 2 

~Vp + (F(pi)V(^i - F{p2)V(t>2)) ■n = Pp and y V(^ + /3<^ 



(3.16a) 
(3.16b) 

(3.17) 
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Using (3.16) we calculate (integrating by parts): 
d 



dt 



n 



(j)p dx = I {(pdtp + pdt(l)) dx 



n 



1 



((F(pi)V0i - F(p2)V02)V</.) dx + ^ I {{E'ipi) - E'{p2))p)dx 
n ^ Jn 

l{F'{pi)\VcPi\^ - F'{p2)\V(t>2\^)p dx 

(F(pi)V0i - F(p2)V02)V,^ - \{F\p{)\V 4>x\^ - F\p2WH^)p\ dx 
n ^ 

+ ll^{E'ipi)-E'ip2))pdx. 



Note that due to the boundary conditions (3.17), all boundary terms that result from the 
integration by parts in the previous calculation vanish. Due to the concavity of F we deduce 
that 

(F(pi)V(/.i - F(p2)V(/>2)V0 - ^(F'(pi)|V(/<i|2 - F'{p2)\V(l)2\^)p = 
= (F(pi) - ^F'(pi)p)|V</)iP + (F(P2) + ^F'(P2)P)|V02|' - (Fipi) + F(p2))V</)iV02 
>^(F(pi) + F(p2))(V(/>i-V(/>2)2. 
Since E is in C^(T) and convex, E' is monotone and we obtain 



d 
dt jn 



dx> I hF{pi) + Fip2)){V(Pi-VcP2f dx. 
In ^ 



We integrate over the interval [0,T] and obtain (since p(x,0) = and (f){x,T) = 0) that: 

/ J(F(pi) + F{p2)){V^i - V02)' dx < 0. (3.18) 

Jn ^ 

The function F is positive, therefore (3.18) implies that 

V(^i = V02 in {pi > 0} U {p2 > 0}. 

2 

Then pi is a solution of dtp = \^p — V • {F{p)V(j)2) = and we conclude (by uniqueness) 



that pi = P2- 



U 



4. Numerical simulations 

In this section we present a steepest descent approach for the solution of the parabolic optimal 

control problem (2.12). This method can be used for convex optimization problems and defines 

an iterative scheme to determine the optimal p and v. 

The presented numerical simulations focus on the different behavior of the classical model of 

Hughes (2.21) and the corresponding general mean field modification (2.22). This comparison 

requires the simulation of (2.21 ), which is done using a finite volume scheme for the nonlinear 

conservation law and a fast sweeping method for the Eikonal equation, as presented by Di 

Francesco et al. in |17j . 

The steepest descent scheme for (2.12) can be written as: 
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Scheme 1. Let po = Po{x), (j){x,t) = (pnix), where (pn is the solution of the Eikonal equation 

G(p), 
Fip) 



in the classical Hughes model (2.21 ) with p{x, t) = Pq{x) for all t G [0, T] and v = -pTfrVc/), be 



the given initial data. Then the steepest descent scheme reads as: 
(1) Solve the nonlinear convection diffusion equation 

2 

dtp=^Ap-V-iF{p)v) (4.1a) 

2 2 

( — —Vp + j)-n = f3ponTE and ( — — Vp + j) • n = on Tat, (4.1b) 



forward in time. System (4.1) is solved implicitly in time. The resulting nonlin- 
ear equation A{p, u) = is discretized using a mixed hybrid discontinuous Galerkin 
(MHDG) method, cf. |20j . and solved using Newton's method. 

(2) Calculate the backward evolution of the adjoint variable (/> = 4>{x, t) using the previ- 
ously calculated density p = p{x, t) in 

- dtc^ - y A0 - G'{p)v . V<A = -If'{p)\v\^ - ^E'ip), 

2 2 

— :r'^4' ■ n ~ P4> = on Te and — p^Vc/) • n = on Vn, 

with an implicit in time discretization and a MHDG method for the spatial discretiza- 
tion. 

(3) Update the velocity via v = v — t{F{p)v — G{p)'V(p), where r is a suitably chosen step 
size. 

(4) Go to (1) until convergence of the function (|2.12a). 



Egger and Schoberl presented a MHDG method for linear convection dominated problems, 
which can be adapted for the nonlinear problems considered. We use the following basis 
functions for the discretization, i.e. 

/), 0, V G P°(T) and Vp, V0 G RT^{T), 

where P^ denotes piecewise constant basis function on the interval T and RT^ lowest order 

Raviart-Thomas basis functions. The Newton iteration in (4.1 ) is terminated, if ||^(p, ^^)||l2(q\ < 

10-6. 

Throughout this section we consider the domain il = [—1, 1] with exits located at x = ±1. 

The interval Q is divided into a set of equidistant intervals of length h. Furthermore the 

maximum density Pmax is set to 1. 

4.1. Comparison of the classical and the mean field type Hughes model. We have 



seen in Section 2.5 that the proposed mean field model has a similar structure as the clas- 
sical model by Hughes. Hence we want illustrate the behavior of both models with various 
experiments. In the first example we focus on the infiuence of the boundary conditions as 



discussed in Section 2.6 In the second example we illustrate the basic difference of the mean 



field game approach and the classical model of Hughes. 

4.1.1. Behavior for different values of (3. We consider the time interval t G [0,3]. In the 
classical Hughes model ( |2.21 ) the boundary conditions are set to 

2 

H±l,t) = and (-y Ap + pf{p)V^) ■n = pp, 
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Figure 1. Solution of the classical Hughes model (|2.21|) with initial datum 

1 

3 



po{x) = 3 at times t = 0.1, 0.7, 1.5 for different values of 



a 



with f{p) = 1 — p. For the corresponding mean field model (2.22) the functions p and 



satisfy (2.24) and (2.26) at x = ±1. In both cases we use the same diffusion coefficient, i.e 
a = 0.1. The spatial discretization is set to /i = 5 x 10~^ in (2.21) and h = 10^^ in (2.22), the 
time steps to Af = 10^^ and Af = 10~^ respectively. The different magnitudes of the time 



stepping can be explained by the explicit in time discretization of (2.21 ) and the implicit time 
discretization of (2.22). Figures [I] and p] show the evolution in time of the solutions of (2.21) 
and (2.22) for different values of (3. Although the models have a very similar structure, their 



behavior is different. In the mean field model small congestions at the boundary are visible 
for /3 = 1. Furthermore we do not observe the immediate vacuum formation at x = as in 



Hughes model (2.21). People rather tend to "wait for a little while" at the center and then 



start to move at a higher speed. The expected equilibration of </> in (2.22) is clearly visible in 
Figure [2] for all values of /3. 

4.1.2. Fast exit of several groups. In this example we consider three groups, which want to 
leave the domain as fast as possible. The particular initial datum is given by 



Po{x) 



A similar example was already considered in [17], where the simulations showed that the a 
small part of the group located around x = 0.6 initially splits to move to the exit at x = — 1, 



0.8 


if -0.8<x< -0.6 


0.6 


if - 0.3 < X < 0.3 


0.95 


if 0.4<x < 0.8 





otherwise. 
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Figure 2. Solution of the mean field type Hughes model (2.22) with initial 



datum po{x) = i at times t = 0.1, 0.7, 1.5 for different values of a 



but later on turns around to take the closer exit at x = 1. This behavior can be explained by 
the fact that the pedestrians in Hughes model ( |2.21 ), adapt their velocity in every time step 
(depending on the overall density at that time). Due to the initially high density of people 
located in front of exit x = 1, parts of the group start to move towards x = — 1, but turn 
around when the density is decreasing as more and more people exit. 



We do not expect to observe this behavior for the modified Hughes model (2.22). One of the 



underlying features of the proposed optimal control approach, is the fact that each pedestrian 
knows the distribution of all other people at all times. Therefore he/she is anticipating the 
behavior of the group in the future and will in this case rather wait than move to the more 
distant exit. This expected behavior can be observed in Figure [3j Here the mean field model 
(2.22) was solved using the steepest descent approach detailed before. The parameters were 



set to: 



10"^ At 



5 X lO^^T 



3 and /3 = 1. 



Note that the optimal control formulation of (2.22) is not a convex problem (the function 



H{p) = p{l — p)^ is only concave for p G [0, |]), but the steepest descent approach converged 
without any problems. For Hughes model we chose the following parameters: 



h = 5x 10"^ At = 5 X 10^^ and /3 = 1. 
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(a) Solution of the classical Hughes model (2.21 \ (b) Solution of the mean field optimal control ap- 
proach (12. 221) 



Figure 3. Fast exit scenario for three groups 



4.2. Linear vs. nonlinear E. In the final example we ihustrate the behavior for different 
functions E = E(p). In particular we set 



F{p) = G{p) 



and E{p) 



3p 

p3p ■ 



The later choice of E actively penalizes regions of high density and urges the group to get 
out more efficiently. In this example the group is initially located at 



Pq{x) 



if -0.25< x<0.4 

otherwise. 



The parameters are given hy h = 10^^, At = 0.1, T = 3 and /3 = 1. Figure [4| shows the 
evolution of the density p for the different choices of E in time. As expected we observe that 
the group exits much faster and spreads out more evenly for E = e^^ than in the linear case. 

Conclusion 

In this paper we presented a novel mean field game approach for evacuation and fast exit 
situations in crowd motion. We motivated the model on the microscopic level and discussed 
its generalization in the mean field limit. The optimality system of the resulting parabolic 
optimal control problem, is a mean field game and establishes interesting links to well known 
models, like Hughes model for pedestrian flow. Furthermore the proposed model gives new 
insights into the mathematical modeling of boundary conditions for pedestrian crowds. We 
present first existence and uniqueness results for the proposed optimal control approach and 
illustrate the behavior of solutions with various numerical simulations. 

The general formulation of the proposed model poses interesting challenges for future re- 
search. In general the derivation of mean-field limit with nonlinear mobilities is of significant 
importance and raised a lot of interest in the scientific community, cf. e.g. [13]. In addition 
challenging questions and problems arise in the mathematical analysis of the optimal control 
approach for convex as well as non-convex problems. Another important direction of research 



MEAN FIELD GAMES WITH NONLINEAR MOBILITIES IN PEDESTRIAN DYNAMICS 



19 



Time t=0.1 



Time t^O.4 




-0.5 0.5 

Time t-0.7 



0.5 0.5 

Time t-1 




Figure 4. Evolution of the density p for different energies in time. The blue 
line corresponds to E = 3p, the red one to E = e^^ . 



will focus on more realistic modeling assumptions. We have seen in the second example of 
Section [4] that people anticipate the density of all others in time, which is questionable in 
many situations. Therefore it would be interesting to study the behavior of (2.12) with a 



temporal discount factor, as proposed in [30] . Finally the development of a numerical solver 
for 2D problems and its generalization to non-convex problems would allow for more realistic 
simulations. 
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